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Abstract 

We calculate the Yang-Lee zeros for characteristic temperatures of the 
helix-coil transition in a continuum model of polyalanine. The distribution of 
these zeros differs from predictions of the Zimm-Bragg theory and supports 
recent claims that polyalanine exhibits a true phase transition. New estimates 
for critical exponents are presented and the relation of our results to the Lee- 
Yang theorem is discussed. 
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1 Introduction 

Recently, there has been an increased interest in the statistical physics of biolog- 
ical macromolecules. For instance, the conditions under which a-helices, a common 
structure in proteins, are formed or disolved, have been extensively studied 0. Tra- 
ditionary, the characteristics of the observed sharp transition between random coil 
state and (ordered) helix, have been described in the framework of Zimm-Bragg-type 
theories ||. In these theories, the molecules are approximated by a one-dimensional 
Ising model with the residues as "spins" taking values "helix" or "coil" , and solely 
local interactions. Hence, thermodynamic phase transitions are not possible in these 
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theories. However, in previous work || f|, [I| evidence was presented that polyala- 
nine exhibits a phase transition between the ordered helical state and the disordered 
random-coil state when interactions between all atoms in the molecule are taken into 
account. Investigating the finite size scaling of various quantities such as specific 
heat, susceptibility and the Fisher zeros, we studied the nature of this phase transi- 
tion and presented estimates for the critical exponents u, a and 7. Here, we continue 
our previous investigation by exploring the helix-coil transition in polyalanine from a 
different point of view. We intend to demonstrate in greater detail the differences in 
the critical behaviour of our model and the Zimm-Bragg theory, and we will present 
additional critical exponents. 

Yang and Lee have established long ago || that the statistical theory of phase 
transitions can also be described by means of the distribution of the zeros of the 
grand partition function in the complex fugacity plane. In the case of the Ising 
model this description corresponds to study this model with an external magnetic 
field H mathematically extended from real to complex values. Although complex 
values for the fugacity are not physical, the characteristics of the phase transition 
can be extracted from the distribution of the zeros. As the volume of the given (fi- 
nite) system increases, the number of complex zeros grows and they will close onto 
the positive real fugacity axis. In the thermodynamic limit, those zeros circum- 
scribe closed regions on this plane, thus defining thermodynamic phases which are 
themselves free of zeros. In Ref. ]7| Lee and Yang proved the famous circle theorem, 
namely that in the Ising model with ferromagnetic couplings the complex zeros lie 
on the unit circle in the complex y— plane, y = e~ 2 ^ H . 

The goal of our paper is to investigate whether the helix-coil transition in bio- 
logical macromoelecules can also be described in the frame work of the Yang-Lee 
zeros. For polyalanine, the analogous of the magnetization in the Ising model is the 
number of helical residues M. Formally we can introduce a (non physical) external 
field H as the conjugate variable to this order parameter M (which plays the role 
of a magnetization), and study the Yang-Lee zeros. Such an investigation allows us 
not only to check our previous results with an independent method, it opens also 
the possibility to calculate new critical exponents and to point out much clearer the 
differences between our all-atom model of polyalanine and the Zimm-Bragg theory. 
If the helix-coil transition can indeed be described by the Zimm-Bragg model, then 
the distribution of the Yang-Lee zeros in function of the temperature should re- 



different distribution would demonstrate that the helix-coil transition in polyalanine 
is not accurately described by the Zimm-Bragg theory. 

Before proceeding, we give the outline of the paper. In the next section we de- 
scribe the numerical evaluation of the partition function. Our approach and numer- 
ical results are presented in the third Section, which is followed by our conclusions. 



2 Methods 

Our investigation of the helix-coil transition for polyalanine is based on a detailed, 
all-atom representation of that homopolymer. Since one can avoid the complications 
of electrostatic and hydrogen-bond interactions of side chains with the solvent for 
alanine (a nonpolar amino acid), explicit solvent molecules were neglected. The 
interaction between the atoms was described by a standard force field, ECEPP/2,|| 
(as implemented in the KONF90 program ||) and is given by: 

E t ot = E c + E LJ + E H b + E tor , (1) 
Ec = (2) 




Eu = £ Rf-5r ■ (3) 



Ehb = E(^-ffe). (4) 
E tor = J2 U i( 1±C0S ( n iXi))- (5) 

I 

Here, rjj (in A) is the distance between the atoms i and j, and xi is the /-th torsion 
angle. Note that with the electrostatic energy term Ec our model contains a long 
range interaction neglected in the Zimm-Bragg theory 0. It was conjectured that it 
is this long range interaction in our model and the fact that it is not one-dimensional 
which allows the existence of the observed phase transition in Refs. ||, |], [5|. We 
remark that the ID Ising model with long-range interactions also exhibits a phase 
transition at finite T if the interactions decay like l/r°" with 1 < a < 2 . 



Simulations of detailed models of biological macromolecules are notoriously dif- 
ficult. This is because the various competing interactions within the polymer lead 
to an energy landscape characterized by a multitude of local minima. Hence, in 

ftip lniff-tpmnprflturp rpcrinn rannniral TVTnntp ("Inrln nr mnlprnlnr Himamins simnln- 



thermalize within the available CPU time. Only recently, with the introduction of 
new and sophisticated algorithms such as multicanonical sampling [T]J and other 



generalized- ensemble techniques ]T2|] was it possible to alleviate this problem in pro- 



tein simulations [O]. For polyalanine, both the failure of standard Monte Carlo 
techniques and the superior performance of the multicanonical algorithm are exten- 
sively documented in earlier work |14|. For this reason, we use again this simulation 



technique for our project where we considered polyalanine chains of up to iV = 30 
monomers. 

In the multicanonical algorithm |ll] conformations with energy E are assigned 
a weight w mu (E) oc l/n(E). Here, n(E) is the density of states. A simulation with 
this weight will lead to a uniform distribution of energy: 

P mu {E) oc n{E) w mu (E) = const . (6) 

This is because the simulation generates a ID random walk in the energy space, 
allowing itself to escape from any local minimum. Since a large range of energies are 
sampled, one can use the reweighting techniques [[15] to calculate thermodynamic 
quantities over a wide range of temperatures T by 



dx A{x) w-\E(x)) e~ pE W 

<A> T = J — r , (7) 

dx w~\E(x)) e-? E{x) 



where x stands for configurations. 

Note that unlike in the case of canonical simulations the weights are not a 
priori known. We needed between 40,000 sweeps (N = 10) and 500,000 sweeps 
(N = 30) for the weight factor calculations by the iterative procedure described in 
Refs. |ni |1(J. All thermodynamic quantities were estimated from one production 
run of N sw Monte Carlo sweeps starting from a random initial conformation, i.e. 
without introducing any bias. We chose N sw =2 000 000, 2 000 000, 4 000 000, and 
3000 000 sweeps for N = 10, 15, 20, and 30, respectively. 

It follows from Eq. (^) that the multicanonical algorithm allows us to calculate 
estimates for the spectral density: 

n(E) = P mu (E)wrJ a (E) . (8) 



We can therefore construct the corresponding partition function for our all-atom 



M at a given energy, 

N-2 

Z(u,y) = £ E^ M W, (9) 

M=0 _E 

where w = e _/3 with /3 the inverse temperature, (3 = 1/ksT and y = e~ H . We have 
absorbed the factor f3 in the definition of the field H which is the conjugate variable 
to our order parameter M. The later quantity is defined by the condition that for 
a helical residue the dihedral angles (0, ip) fall in the range (—70 ± 20°, —37 ± 20°). 
Note, that for convenience we have discretized the energy in Eq. (|9|). In previous 
work we found that the partition function zeros depend only weakly on the energy 
bin size, and for the present study we choose energy bins of length 0.5 kcal/mol. 



3 Yang-Lee Zeros of Polyalanine 

In order to analyse the helix-coil transition in our all-atom model of polyalanine we 
first determine the (pseudo) critical temperatures. This can be done in various ways. 
For instance, in Ref. Q these temperatures where calculated from the distribution 
of the Fisher zeros [ I7| . Here, we check our previous results by calculating estimates 
for the pseudo-critical temperatures in a different way. 

Our temperature-driven transition model can be analysed by considering the 
coexistence of two phases: a disordered (coil) phase and an ordered (helix) phase at 
the critical temperature T c . This behaviour can be described by the existence of an 
equilibrium probability distribution P^{M) for the corresponding order parameter 
M []18 , l9| . This approach leads to a definition of a finite size critical temperature 



by approximating the equilibrium distribution by a sum of two distributions, each 
one characterizing the corresponding phase. This ansatz gives origin to a two-peak 
like distribution separated by a minimum and has led to the finite size scaling 
study of the coexistence of bulk phases P0| . Following these works, we consider 



for our purpose the histograms of the helicity for all chain sizes N, which is our 
order parameter (and the equivalent to the magnetization in the Ising model), at 
different temperatures j3. It is straightforward to calculate the histograms of helicity 
distribution at these temperatures from our multicanonical estimates of the density 
of states, 

W (M) = £n(£,M)u B ( 10 ) 



predominant single peak at low (high) values of the helicity. Hence, we define 
our critical temperatures by the condition that the reweighted histogram of the 
magnetization have two peaks of equal heights. As an example we show in Fig. 1 
the histogram normalized to area one, for N = 30 at T c = 518 K. 

Table 1 summarizes the so obtained set of critical temperatures by this method. 
The quoted values and errors were obtained by reweighting independent binned 
data to temperatures T c *, i = 1, 2, nbins, where one obtains equal heights in 
the histograms for our order parameter M. Our statistics were based on nbins = 
19, 20, 16 and 3 for N = 10, 15, 20 and 30, respectively, for comparison we also list our 
previous estimates [|J which were derived from an analysis of the Fisher zeros and 
the position of the maximum of specific heat C v . We note, that our new results agree 
within the errors with our previous estimates for the pseudo-critical temperatures 
which demonstrates the validity of our approach. 

After having demonstrated the reliability of our data we calculate now at the so 
obtained critical temperatures the distribution of the Yang-Lee zeros by means of 
Eq. (|9|). Our results are displayed in Fig. 2 for the case of chain lengths iV = 30. 
In that figure we draw in addition a unit circle to point out the relative position 
of our zeros around it. Similar distributions are obtained for smaller systems. The 
presented error bars for the zeros in this figure, were estimated by considering the 
new partition functions at the limiting temperatures T C (N) ± AT C (N) where AT C (N) 
is the error in our pseudo-critical temperatures T C (N). Our error bars clearly exclude 
the possibility that the complex zeros are distributed on the unit circle. On the other 
hand, the Yang-Lee zeros for the ID Ising model (which can be calculated exactly) 
fall on the unit circle. This is shown in Fig. 3 for chain lengths N = 5, 15 and 30.0 

After having demonstrated the difference in the two models at the critical point, 
we extend our study to the behaviour of our polyalanine zeros as function of the 
temperature. In Fig. 4 we show for our smallest polyalanine chain how the zeros 
change with temperature u(f3), measured in multiple of the critical value u c . Again, 
we draw for illustration a unit circle to call attention to the relative position of 
our zeros. We exhibit the smallest chain since the density of zeros increases with 
growing chain length making it difficult to the eyes to trace the movement of the 
zeros in the complex y-plane. As the temperature increases, u((3) > u c , we observe 
that the zeros move away from the unit circle and consequently the edge zero, the 
zero closest to the positive real y— axis, moves away from the real axis. A similar 



behaviour is known for the g-state Potts model with q > 2 on finite lattices for 
temperatures larger than the critical one [2^. For u{(3) < u c , we find again that 



the zeros do not fall all on the unit circle. Hence, the Lee- Yang theorem does not 
apply for our model demonstrating again that polyalanine should not be described 
by Ising-type models. 

Figure 4 indicates that the angle 9q, the angle at which the edge zeros deviate 
from the real y-axis, grows as the temperature is increased above T c , i.e. u(f3) > u c . 
For larger lengths N, the number of zeros increases and the edge zeros move towards 
the real axis. Note that the movement for the edge zero seems to follow an unit arc, 
leading to Re(y) = 1 and the real value H = in the thermodynamic limit. The 
movement of the edge zeros towards the real axis itself is a signal for the possibility 
of a true phase transition in the thermodynamic limit if the locus of zeros cut the 
positive real axis at finite temperature. In that case one would expect the following 
scaling relation for the edge zeros ||24||: 



6 (t,N) = N- yh/d 9 (tN yt/d ) , (11) 

where t = (T — T c )/T c , yh is the magnetic scaling exponent and yt = 1/v is the 
thermal scaling exponent. Note, that we have introduced the scaling variable N 1 ^ 
in place of the usual linear length L. This is because the number of monomers N is 
the natural quantity to describe polyalanine chains. Since we have no theoretical in- 
dication to assume a particular integer geometrical dimension d for our polypeptide, 
we present our results as a function of this parameter d. 

Now, assuming that the above finite size scaling relation holds for polyalanine 
chains we can calculate the "magnetic" scaling exponent |T7 



y h /d=(j3 + y)/du (12) 

if we take t = 0, i.e. u{(3) = u c . For this reason we display in Fig. 5 the linear 
regression for In (9 (N)). The least-square fit gives us 

y h /d= 0.97(4), (13) 

with the goodness of fit Q = 0.67. This value for Q means the fit is acceptable for 
the present precision data and corroborates our expectation. 

If we take our estimate in Eq. (O) and the reported values for 7 and dv from 



~r> „f H7TT 1 ne/1 A \ J,. n no/cr\ l j-„ „ a a 1 



precise since it envolves three parameters (j3 = —0.16(15)), while the second one, 
relying on two previous estimates, gives (3 = —0.065(74). As one should expect, both 
estimates agree with each other within the error bars demonstrating the consistency 
of this study with earlier investigations |3|, f|, |5|]. Note also that our estimates for the 
critical exponent (3 do not help us to answer a question we had to leave unanswered 
in the earlier work, namely, whether the helix-coil transition is a first order phase 
transition {(3 = 0) or a (strong) second order phase transition (J3 ^ 0). This is not 
surprising since clarifying the nature of a phase transition can be a challenging issue 
even for simpler systems (such as Potts models |?5|1 , to name only one example) 
and requires dedicated large-scale simulations. The important point here is rather 
that our model of polyalanine, despite its complexity, can be described in terms of 
a critical theory. 

4 Conclusion 

We have extended the Yang-Lee zeros analysis to an all-atom model of polyalanine 
polymers of various chain lengths. Our results were compared with the Zimm- 
Bragg model which is often used to describe helix-coil transitions in Biopolymers. 
Our analyses of the Yang-Lee zeros for both models give clear evidence for our claim 
that the helix-coil transition in an all-atom model of polyalanine is not adequately 
described by the Zimm-Bragg theory. Our results confirm earlier work where we 
pointed out the possibility of a thermodynamic phase transition in polyalanine which 
is forbidden in the Zimm-Bragg model. 
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Table: 



Table 1: Critical temperatures Tj} lst as obtained by our analysis for polyalanine 
chain lengths N. For comparison, we show also previous estimates obtained from 
both an analysis of Fisher zeros (j' Flsher ) and the location of the peak in the specific 
heat (7f«). 
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Figure Captions: 



Figure 1. Histogram for the number of helical residues in conformations for 
chain length N = 30. 

Figure 2. Distribution of complex partition function zeros in the y— plane with 
error bars, at T C (N = 30). We draw an inexistent unit circle to show how the zeros 
are distributed around it. The zeros do not circumscribe an unit circle as expected 
for the Lee- Yang theorem in Ising like systems with positive interactions. 

Figure 3. The Yang-Lee zeros for the Ising chain with lengths N = 5, 15 and 30. 

Figure 4. The behaviour of the partition function zeros for polyalanine with 
chain length N = 10 in function of the temperature u(f3). 



Figure 5. Linear regression for ln(#o), Eq. flTT|), at the critical temperature. 
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Fig. 1: Histogram for the number of helical residues in conformations for chain 
length N = 30. 
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Fig. 2: Distribution of complex partition function zeros in the y— plane with 
error bars, at T C (N = 30). We draw an inexistent unit circle to show how 
the zeros are distributed around it. The zeros do not circumscribe an unit 
circle as expected for the Lee- Yang theorem in Ising like systems with positive 
interactions. 
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. 3: The Yang-Lee zeros for the Ising chain with lengths iV = 5, 15 and 30. 
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4: The behaviour of the partition function zeros for polyalanine with chain 
th N = 10 in function of the temperature u(f3). 
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Fig. 5: Linear regression for ln(0o), Eq. (|TT|) , at the critical temperature. 



